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Off-diagonal Wave Function Monte Carlo Studies of Hubbard Model I 
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We propose a Monte Carlo method, which is a hybrid method of the quantum Monte Carlo method 
and variational Monte Carlo theory, to study the Hubbard model. The theory is based on the off-diagonal 
and the Gutzwiller type correlation factors which are taken into account by a Monte Carlo algorithm. In 
the 4x4 system our method is able to reproduce the exact results obtained by the diagonalization. An 
application is given to investigate the half-filled band case of two-dimensional square lattice. The energy 
is favorably compared with quantum Monte Carlo data. 
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I. Introduction 

Strongly correlated electron systems have been inves- 
tigated by many researchers in order to understand the 
mechanism of superconductivity of the high-T c oxide su- 
perconductors. The effect of strong correlation between 
electrons is important for high-T c superconductivity and 
many unconventional phenomena such as metal-insulator 
transition and heavy fermions with a huge large mass. 
The Hubbard model is a basic model for strongly corre- 
lated electrons in metals. The Hubbard model has been 
regarded as a model to describe the Mott transition in 
compounds such as NiO and MnO. 1 It is considered that 
the Hubbard model contains basic physics of high-T c su- 
perconductivity as a simplified one-band model of the 
three band Cu-0 model. 2 A possibility of superconduc- 
tivity for the 2D Hubbard model has been controversial 
for many years. 3-5 The study of the Hubbard model pro- 
vides us insights having important implications concern- 
ing the origin of high-Tc superconductivity. The possi- 
bility of the superconducting state is recently reported 
for the Hubbard model 6-9 and the two-chain Hubbard 
model. 10-12 The phase diagram is still far from well un- 
derstood in two dimensions. 

The quantum Monte Carlo method is a method to treat 
the strong correlations exactly. However, the applicabil- 
ity is restricted to moderately correlated region because 
of a sign problem. The variational Monte Carlo method 
has a feature characterized by a wide applicability from 
weak to strong correlation. The Gutzwiller wave function 
is a standard trial wave function for itinerant correlated 
electrons. In the Gutzwiller function only the correlation 
at the same site is taken into account and thus the wave 
function is very simple in its form. An analytic evalua- 
tion of expectation values is, however, very difficult for 



the Gutzwiller function in spite of its simplicity. This dif- 
ficulty is overcome by a Monte Carlo method and many 
calculations have been performed using the Monte Carlo 
algorithm. 6,13-15 

It is revealed that the Gutzwiller function will give 
higher energies compared to exact values. Thus we think 
that the Gutzwiller function should be improved to inves- 
tigate the ground state more precisely. The purpose of 
this paper is to propose a hybrid method of the quan- 
tum Monte Carlo (QMC) calculations and variational 
Monte Carlo (VMC) theory. Following an idea of QMC, 
a trial function is improved and the expectation val- 
ues of energy and other quantities are calculated using 
VMC for the parameters optimizing the energy expecta- 
tion value. The off-diagonal wave function Monte Carlo 
method (OWMC) in this paper deal with Gutzwiller pro- 
jection and off-diagonal correlation operators taken into 
account by a Monte Carlo algorithm. We expect to be 
able to extrapolate correct expectation values from the 
data obtained for off-diagonal wave functions. 

We show advantages of OWMC in the following: (1) 
The calculations give an upper bound of the exact en- 
ergy since OWMC is based on the variational theory. (2) 
Variational wave functions are improved systematically. 
(3) There is no sign problem. (4) The large- [/ systems 
are tractable. (5) Introducing the order parameters we 
can investigate ground state properties characterized by 
the long-range ordering. 

The paper is organized as follows. The section 2 is 
devoted to present the Hamiltonian and wave functions. 
The method of calculations is also shown. The following 
two sections are assigned to a description of results. The 
summary and discussion are presented in the final section. 
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II. Hamiltonian and wave functions 
A. Wave functions 

The Hamiltonian is given by the Hubbard model: 

H = -tJ2( c L c 3<? + h - c +U^2n n n il , (1) 

(ij) 1 

where i%i a — c^Cia is the number operator and U is the 
strength of on-site Coulomb interaction. The energy is 
measured in units of t throughout this paper and the 
number of sites is denoted as N. In this paper the Hub- 
bard model is considered in two space dimensions. The 
Coulomb interaction U is expected to bring about the 
metal-insulator transition, antiferromagnetic order and a 
possibility of the anisotropic superconducting state. 
Let us start with the Gutzwiller function given by 

ip G = Pc^a = exp(-a^n 4T n a )^o = ^ (0) , (2) 

i 

where Vo is the free fermion ground state. Pq is the 
Gutzwiller operator given by Pg = — (1 ~ 9) n il n il) 
where g is the variational parameter in the range of < 
g < 1 and a = log(l/g). Our method is based on the 
fact that the ground state eigenfunction is written as in 
QMC 

j> = e-^o * e^ K e-^ v ' ■ ■ ■ e~^ K e~^ v ' ^ 

(3) 

for large j3 = e\ H — ■ + e n and small e, (i = 1, ■ • • , n). H 
is written as H = K + V where K denotes the kinetic 
energy part and V' = UV denotes the on-site interaction 
part: V = J2i n i1 n il- Since the last factor e~ eV ?Ao m 
eq.(3) is regarded as the Gutzwiller function, one way 



to improve ipG is achieved by multiplying e XK to ipc- 
Therefore we can consider 16 

V> (1) = e- XK e- aV iP . (4) 

It has been shown that the Gutzwiller function is im- 
proved appreciably by the off-diagonal correlation factor 
g-A-R: 16-18 rp ne nex ^ s ^ e p j s ^ multiply e~ a v in order to 

control the double occupancy in ip^K We can again op- 
erate e~ XK to improve the wave function. A second-level 
wave function is given by 

^ = e- x ' K e- a ' v e- XK e- aV TPo, (5) 

where A, A', a and a'(= log(l/</)) are variational param- 
eters. Next wave function is written as 

^(3) = e -X"K e - a "V e -X'K e - a 'V e -XK e - a V % (g) 

A, A', A", a, a' and a" are variational parameters. It is 
considered that with optimized variational parame- 
ters approaches the ground state wave function in eq.(3) 
as m grows larger. Although a sign problem will occur for 
large to, the sign never brings about a problem for small 
to. If we can extrapolate the expectation values from the 
data obtained using we can estimate exact 

values within Monte Carlo errors. 

B. Method of calculations 

A Monte Carlo algorithm developed in the auxiliary- 
field quantum Monte Carlo calculations 19 enables us to 
evaluate the expectation values for the wave functions 
tp^ , ip( 2 \ ■ ■ ■. Using the discrete Hubbard-Stratonovich 
transformation, the Gutzwiller factor is written as 



exp(-a^n iT n a ) = (1/2)^ ^ exp[2o ^ S;(n iT - n tL ) - -^(n 4 T (7) 

i Si i i 

where a is given by cosh(2a) = e"/ 2 . N denotes the number of sites. The auxiliary fields Sj takes the values of ±1. 
The norm (V'gIV'g) is written as 

WgWg) = (V2) 2Ar J2 nW|cxp(^( U ))exp(^( S ))|^}, (8) 

where the potential h a {u) is given by 

h a (u) = 2aa S ^u l n la - ^ ^n i(T . (9) 

i i 

Then the weight is written as a sum of determinants, 20,21 

(^ G |Vg)= const. (l/2) 2N Yl n det ^o t exp(U CT ( U , a ))exp(y ff ( S ,«))^), (10) 

{«,}{««} ° 

where V a (s, a) is a diagonal N x N matrix corresponding to h° (s) written as V{s, a) = diag(2acrsi — a/2, • • • , 2aasiy — 
a/2), where diag(a, • • •) denotes a diagonal matrix with elements given as the arguments a, • • •. For the free fermion 
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state, the elements of <py are given by plane waves: (</>q )ij = exp(ii"j • kj) (i = 1, • • • , N; j = 1, • • • , N e /2 where N e is 
the number of electrons). We may take k, below the Fermi surface corresponding to the free fermion ground state as 
a trial state. To represent 4>o as a re& l matrix, we take {4>Z)ij as cos(rj • kj) or sin(rj ■ kj). A more general choice is 
possible by incorporating the spin-density wave (SDW) long range order. The one-particle antiferromagnctic ordered 
state is given by 



V'sdw - II( u k c k T + W k c k+Q T ) Il^k'ii - ^k'i+QP^' 



(11) 



where u k = [(l-e^/(e^ + A 2 AF ) 1 / 2 )/2] 1 / 2 and u k = [(l + e k /(e k + A 2 ^) 1 / 2 )^] 1 / 2 . A AF is the antiferromagnetic order 
parameter optimizing the energy. Q denotes SDW wave vector given as Q = (n, ir). The elements of </>q corresponding 
to Vsdw are written as 



(4>%)ij = u k .exp(iri • kj) + signcrv k exp(ir l • (kj + Q)). 



(12) 



In real representation they are given by u k cos(ivkj) +sign<rw k cos(rj • (kj+Q)) and u k sin(rj-kj) +sign<7z; k sin(r, 
(kj+Q)). 

Similarly the norms including the off-diagonal factors are written as 



= const.(l/2) 2JV ^ f[ det^cxp^K a))exp(-\K a )exp(-\K a )exp(V a (s, a))^) 



(13) 



and 



(V> (2 V (2) ) = const. (1/2) 4JV Y[det(^exp(V a (u,a))cxp(-XK a )cxp{V a (v,a')) 

{u 3 }{v e }{t k }{ Si } <? 

x exp{-\' K a )cxp{-\' K^^xp^ 17 (t, a'))exp(-XK' 7 )exp(V a (s, a))$J), 



(14) 



where K a is a matrix corresponding to the kinetic part 
of the Hamiltonian: 



(K a )ij = —t if are nearest neighbor sites, 



In order to evaluate the expectation value we generate 
the Monte Carlo samples by the importance sampling 
with the weight function = where 



= otherise. 



(15) 



w„ 



For two states \ip a ) and \il} a ) represented by Slater de- 
terminants <jf and <p a as defined above respectively, the 
single-particle Green function is written as 21 



det(</>^exp(F ff (u, a)) ■ ■ ■ exp(V(s, a))<j>%). 



(17) 



(16) 



for (V'lV'' 7 ) + 0. 



When we update the Ising variable from the old Sj to the 
new one s\, the ratio of is calculated to determine 

whether to accept or reject the new configuration. In 
the process of updating Sj to s\, the Gutzwiller factor 
exp(V CT (s, a)) is written as 



exp(V" CT (si 



1 1 ^27 



, s N 



,a)) = diag(e 2aCTSl - a / 2 ,---,e 2aCTS '- Q/2 ,---,e 2a,TSN - Q/2 ) 



= (1 + A tT )exp(V' T (s 1 ,---,s i ,---,s N ,a)), (18) 
where A a is a diagonal matrix given by A a = diag(0 • • • 0, e 2aCT ( s ;- s i) — 1 7 • • •). Then the ratio of \w a \ is given by 20 

r a = |det(<Ao t exp(y <7 (u, a)) ■ ■ ■ (1 + A a )cxp(V a (s, a))^)|/|det(^ t exp(y CT ( U , a)) • • • exp(V°(s, a ))</%)\ 
= |det(L(l + A ff )i?|/|det(Li?)| 

= |det(l + LA a RJ)\ = |1 + {A a )ii(G rcr )ii\, (19) 
where the right Green function is written as 

G ra = cxp(V a (s,a))<j>°jfj cxp(V a (u,a))exp(-\K a ) ■ ■ ■ exp(-\K a ), (20) 
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with 

J = {(f>' T Jexp(V a (u,a))exp(-XK' T ) ■ ■ ■ exp(-XK a )exp(V a (s, a))#f} -1 . (21) 



i 

The updated Green function is straightforwardly calcu- We continue the updating procedure for other Ising vari- 
lated from the old one if the new variable s- is accepted: ables. The relation in eq.(22) is derived as follows. G^ w 

is written as 

GZ W = (1 + A a )[G ra - G r „A a (l + G rcr A a )~ 1 G riJ }. 

(22) 



GZ W = (l + A a )exp(V^s,a))^J'^exp(V (r (u,a))---exp(-XK^ = (I + A a )G™ w , (23) 
where J' is given by 

J' = {$}exp{V a {u, a))exp(-XK°) ■ ■ ■ cxp(-A^)exp(T/ CT ( Sl , ■ • • , «<, • • • , s N , a))^}" 1 . (24) 
Inserting the relation 

1 = jfJcxp(V' T {u,a))cxp(-XK< 7 )---cxp(-XK< T )cxp(V' T (s,a))<P° 

= J[(J') _1 - <p a J expiV 7 (u, a))cxp(-XK a ) ■ ■ ■ cxp{-XK <J )A a cxp(V a (s, a))(j>%], (25) 

into the left of J' in cq.(23), G™ w is given as 

Then the eq.(22) is followed. When we update the variable Ui in the left part of w a to new one w-, the updated Green 
function is similarly calculated as 

GZ W = [Ge„ - G ta {l + AtjGitjY 1 A a Gtcr](l + A„), (27) 
where the old Green function is given as 

G la = cxpi-XK 17 ) ■ ■ • cxp(-A J ^ <T )cxp(F ,T (s, a))<P°J<f>^ cxp(V a (u, a)), (28) 

and A a is a diagonal matrix = diag(0, • • • , 0, e 2a <?«-^) - 1, 0, • • •)• 

i 



Since the Monte Carlo samplings are generated with the term is given by 
weight |tu|, the expectation value is calculated from the 

summation with the sign of w (denoted as signw) . For ex- < c >= N ^(G ) ••signio/ N sigm/; (29) 

ample, the expectation value of nearest neighbor transfer m ~ 

where 



(Ga)ji = (exp(-X"K°) • • • exp(-A^ CT )exp(T/ CT ( S , a)W 9^ ' exp(V° (u, a)) 

x cxp(-\K a ) ■ ■ • exp(-X"K a ))ji. (30) 



J2 m denotes a summation over the Monte Carlo samples. III. Comparison with exact results and other QMC 

The two body correlation functions are similarly calcu- methods 
lated using the Wick's theorem. 

In standard projector QMC approaches, we encounter 

an inevitable sign problem. In our calculations it has This section is devoted to examine the validity of our 

turned out that the negative sign problem is not serious, method by comparing energies obtained using the off- 

which enables us to estimate the expectation values even diagonal functions with the exact results for the 2D Hub- 

in two space dimensions. bard model. In Table I we compare our results with 
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TABLE I. Ground state energies for the 2D Hubbard model. The boundary conditions are periodic in both directions, ipo is 
set to be the normal free electron state. The Constrained Path Monte Carlo (CPMC) results are from Ref.21. The projector 
Quantum Monte Carlo (QMC) results are from Ref.22. 







U 
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Extrapolated 


CPMC 


QMC 


Exact 


4x4 


10 


4 


-19.39 


-19.57 






-19.582 


-19.58 


-19.581 


4x4 


10 


8 


-17.06 


-17.47 


-17.49 


-17.51 


-17.48 




-17.510 


10x10 


82 


4 


-106.7 


-109.1 


-109.3 


-109.6 


-109.55 


-109.7 





-17.0 



-17.2 

E 

-1 7.4 




-17.6 1 1 

0.4 0.8 

1/(m + 1 ) 

FIG. 1. The energy for N e = 10 and U = 8 on the 4 x 4 
lattice. The results for tpQ, ip^ and ip^ are shown. 

The diamond denotes the exact value. 



TABLE II. Correlation functions for the 4x4 Hubbard 
model with periodic boundary conditions in both directions. 
N e = 10, U — 4 and V>o is the normal free electron state. 







^) 


CPMC 


Exact 




0.808 


0.73 


0.729 


0.733 


C(7T,7r) 


0.460 


0.52 


0.508 


0.506 


A w (l) 


0.080 


0.076 




0.077 


A yy (2) 


0.007 


0.006 




0.006 


A,ir(0) 


0.133 


0.12 




0.122 


A*„(l) 


-0.013 


-0.015 




-0.014 


s(0,0) 


0.525 


0.53 




0.533 


«(1,0) 


-0.098 


-0.091 




-0.0911 


c(0,0) 


0.33 


0.33 




0.326 


c(l,0) 


-0.047 


-0.056 




-0.0539 



the exact values and available data from quantum Monte 
Carlo method for the closed shell case. 

Obviously i/>( m ) (m = 2,3) are improved appreciably 
compared with the Gutzwiller function showing energies 
which are very close to the exact values. Extrapolated 
energies are obtained through an extrapolation from the 
data obtained for ipW = ipQ, ipW , -0' 2 ' and ip( 3 \ as 
is shown in Figs.l and 2. In the above calculations for 
the closed shell case, the lower level wave function ip^> 
is found to be a nice approximation for the ground state 
function. In Table II we show results for the correlation 
functions comparing them with the exact values. 



-107.0 

E 

-108.0 
-109.0 




-110.0 i 1 1 1 1 

0.4 0.8 

1/(m + 1) 

FIG. 2. The energy for N e = 82 and U = 4 on the 10 x 10 
lattice. The results for tpc i}> , ^ 2%> and ip^ are shown. 
The square denotes the expectation value by CPMC and the 
triangle denotes that by QMC. 

In Table II the correlation functions are defined by 

s(i x ,i y ) = (("ot _ n oi)(mi ~ nn)) for r i = (*a,*y), 

(31) 

c(i x ,i y ) = ((n - (n ))(n t - (n,))) for = (i x ,i y ), 

(32) 

where n% — + n^. S(q) and C(q) are Fourier trans- 
forms of s(i x ,i y ) and c(i x ,i y ), respectively: 

S(q) =^exp(iq-r i )s(r i ), (33) 

i 

C(q) = ^exp(iq-r i )c(r i ). (34) 

i 

The superconducting correlation function A a p(£) is de- 
fined by 

A a p(i) = (A{(i + lx)A p (i)), (35) 

where A a (i) = Ciic i+a i - c^Ci +a i (a = x or y; x and y 
indicate an unit vector in x and y direction, respectively). 
The expectation values of correlation functions agree with 
the exact values considerably well. In these calculations 
ip^ satisfactorily reproduces the correlation functions as 
well as the ground state energy. In Table II, parameters 
are given by g = 0.30, A = 0.07, g' = 0.36 and A' = 0.13 
(for U = 4). 
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-14.6 - 



-15.0 - 

E 

-15.4 - 




-15.8* ' ' ' ' 

0.0 0.4 0.8 1.2 

1/(m+1) 

FIG. 3. The energy as a function of l/(m + 1) for N e = 14 
and U = 4 on 4 x 4. For the upper and lower curves an 
initial function ipo is set to be the Fermi sea or SDW state, 
respectively. The diamond indicates the exact value. 



-9.40 - 



-9.60 - 



E 




-10.2 i 1 1 1 1 

0.2 0.4 0.6 

1/(m + 1 ) 

FIG. 4. The energy as a function of l/(m + 1) for iV e = 14 
and t/=12on4x4. For the upper and lower curves an 
initial function tpo is set to be the Fermi sea or SDW state, 
respectively. The diamond indicates the exact value and the 
triangle indicates CPMC result. 



Now let us show an another example where we con- 
sider the 4x4 lattice with 14 electrons for the periodic 
boundary condition in both directions. A famous sign 
problem occurs in this case for which the standard pro- 
jector QMC method gives poor estimates of energy for 
large U. The quantum number of total momentum for 
i/jq can take (0, 7r) and (tt,0) as well as (0,0). Two kinds 
of initial trial states ipQ are examined; one is the Fermi 
sea with zero total momentum for both spin states and 
the other one has non-zero total momentum (0, it) for 
each spin state. We incorporate the antiferromagnetic 
order into a trial wave function in the latter case since 
the energy gain is appreciable. In Table III we show our 
data. The extrapolated values are obtained as shown in 
Figs. 3 and 4 where the solid curve and dotted curve in- 
dicate the energies obtained for SDW and normal initial 
states, respectively. Our calculations are supported by 
the property that the extrapolated values for different 
initial states coincides with each other within statistical 
errors. The energy obtained for ijj^ is comparable with 
that by CPMC for large U and extrapolated values are 
fairly close to exact values. 

IV. Results for the half-filled case 

In this section the 2D square lattice with the half- 
filled band is investigated as a first step toward more 
general cases. The half- filled band systems have been 
studied numerically by the QMC simulation 24 and VMC 
calculations. 16 ' 25 In our calculations we impose the peri- 
odic boundary condition in one direction and antiperiodic 
condition in the other direction to have a unique ground 
state. Two sizes, 6x6 and 8x8, are investigated in this 
section. 

The energy expectation values estimated by the off- 
diagonal wave function Monte Carlo method (OWMC) 
are shown as a function of l/(m + 1) for 6 x 6 in Fig. 5. 
The solid curves are drawn by the least squares method 
for three data points obtained by tp^ , ip( 2 ^ and ip^ . The 



extrapolated values are used to obtain the total energy as 
a function of U as shown in Fig. 6 for 6x6. In Fig. 6 the 
energies for i/jq = Pcipo, Pg^Psdw, ip^ and extrapolated 
values are shown together. For the Gutzwillcr function 
the energy is greatly lowered when the antiferromagnetic 
order is introduced. We obtain lower energies than those 
for -Pg^Sdw by using the off-diagonal wave functions. In 
order to find optimum parameters we use the correlated 
measurements method 26,27 to find the most descendent 
direction in the parameter space. We shift parameters 
along that direction starting from a set of parameters. We 
show our parameters in Table IV. The extrapolations for 
8 x 8 as a function of l/(m + 1) are shown in Fig. 7 where 
the upper and lower curves correspond to U = 8 and U — 
4, respectively. The energy expectation values for 8x8 are 
also shown in Fig. 8 where the energies by ipG, Pg V'sdw, 
ip( 2 ' and extrapolation are compared to one another. The 
variational energy of i/^ 3 ' is very close to that of tp^ 
and the extrapolated value. The results by the quantum 
Monte Carlo simulation by Hirsch are also shown as a 
reference. Although the QMC gives slightly higher energy 
for U — 8, a good agreement between two calculations 
support our method as well as the QMC simulation. It 
is fortunate that the total energy for the intermediate 
value of U is almost independent of the system size, which 
means that the extrapolated energies per site for 6x6 
and 8x8 coincide with each other. Thus the energies in 
Figs. 6 and 8 are expected to agree with the behavior of 
infinite systems. 

The spin structure factor defined as 

jt 

Tin)), ( 36 ) 

can also be evaluated with our method, where N denotes 
the total number of sites. We show S(q) calculated by 
ip^ and ip^ for U = 8 in Figs. 9 and 10, respectively, 
as a function of q x and q v . The strong antiferromagnetic 
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TABLE III. Ground state energies for the 4x4 Hubbard with N e = 14 for the periodic boundary conditions in both directions, 
where (7 = 4 and U = 12. ipo is the non-interacting Fermi sea (normal) or SDW state (SDW) as indicated. The exact 
diagonalization results are from Ref.23. 



u 




4>w 




Extrapolated 


C PMC 


Exact 


4 


normal -14.435 


-15.29 


-15.40 


-15.73 


-15.73 


-15.74 


4 


SDW -15.34 


-15.65 


-15.67 


-15.73 




-15.74 


12 


normal -7.80 


-9.50 


-9.58 


-9.98 


-9.696 


-10.048 


12 


SDW 


-9.75 


-9.80 


-10.0 




-10.048 


TABLE IV. Variational parameters used in the simulation for 6x6. 




U wave function 
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g' 


A' 


g" 


A" 




4 


0.55 
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1.0 


0.0 


1.0 


0.0 






0.18 


0.118 


1.0 


0.0 


1.0 


0.0 




^(2) 


0.47 


0.00064 


0.44 


0.123 


1.0 


0.0 




^(3) 


0.40 


0.105 


0.38 


0.126 


0.398 


0.128 
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0.24 


0.0 


1.0 


0.0 


1.0 


0.0 




^(1) 


0.064 


0.103 


1.0 


0.0 


1.0 


0.0 




^(2) 


0.17 


0.00247 


0.14 


0.103 


1.0 


0.0 




^(3) 


0.18 


0.159 


0.128 


0.197 


0.142 


0.096 



-0.8 



1 1 


1 1 


i 




, -• ' 

U = 8 
U = 4 




. « • • " " 

1 1 1 1 1 



0.0 0.4 0.8 1.2 

1/(m + 1 ) 

FIG. 5. The energy as a function of l/(m + 1) for 6 x 6 at 
half-filling. The upper and lower curves correspond to U = 8 
and (7 = 4, respectively. 




-2.0 ' 1 1 1 1 1 

0.0 4.0 8.0 12.0 

U 

FIG. 6. The energy as a function of U for 6 x 6 at half-filling. 
From the top the energies for tpa, -PgV'sdw, V'' 2 ' an d extrapo- 
lated values are shown. The QMC results for 8x8 are shown 
by open circles as a reference. 



-0.4 



-0.6 - 

lD 

-0.8 - 



-1.0 1 1 

0.0 0.2 0.4 0.6 

1/(m + 1) 

FIG. 7. The energy as a function of l/(m + 1) for 8 x 8 at 
half-filling. The upper and lower curves correspond to U = 8 
and (7 = 4, respectively. 

feature is shown in the figures with highly enhanced com- 
ponent at q = (tt, 7r). 

The sublattice magnetization m has been calculated as 
a function of U with the optimum variational parameters. 
The sublattice magnetization is defined as 

^=^lD- 1 ) J '^T-n i i)| (37) 

3 

m is plotted in Fig. 11 as a function of U where the cal- 
culated values obtained by (for 6x6 and 8x8), 
-PgV'af 25 are shown comparing with QMC method. 24 
Our results are substantially lower than the VMC predic- 
tions with the antiferromagnetically ordered state Pg^af 
and show a good agreement with QMC results. The spin 
fluctuations appreciably reduce the sublattice magnetiza- 
tion but do not destroy the order in two dimensions. 
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FIG. 8. The energy as a function of U for 8 x 8 at half-filling. 
From the top the energies for tpa, PcPsnw, 4> an d extrap- 
olated values are shown. The QMC results are shown by open 
circles. 




FIG. 9. The spin structure factor for ip (2) . U = 8 for 8 x 8 
where q(x) — q x /ix and q(y) = q y /n. 



V. Summary 

We have presented a new variational quantum Monte 
Carlo method for the Hubbard model. Our method is 
based on the property that we can systematically improve 
a variational wave function starting from the Gutzwiller 
wave function. It is remarkable that the variational en- 
ergy is greatly lowered by the Gutzwiller and off-diagonal 
projection operators. For the 4x4 lattice our results 
agree remarkably well with the exact values obtained by 
the exact diagonalization. Our calculations include the 
case where a sign problem occurs for the standard pro- 
jector QMC. The results for 4 x 4 and 10 x 10 (in Table I) 
suggest that lower level off-diagonal functions show up as 
good a level of approximations as the constrained path 




FIG. 10. The spin structure factor for ip (a) . U = 8 for 8 x 8 
where q(x) = q x /n and q(y) = q y /n. 

QMC. The extrapolated values are favorably compared 
with exact results. For the 2D half-filled band system, 
we have calculated the exact energy as a function of U, 
which is lower than that by the antiferromagnetically or- 
dered Gutzwiller function and mostly independent of the 
system size. There is a good agreement with the results 
by the QMC simulation concerning the energy and spin 
structure factor S(q). Thus we have demonstrated that 
the off-diagonal wave function Monte Carlo (OWMC) is 
useful to investigate the ground state for strongly corre- 
lated electrons. 

In our simulation we spent most of time to search varia- 
tional parameters optimizing the energy expectation val- 
ues. We use the correlated measurements method to find 
a most descendent direction along which we shift the vari- 
ational parameters starting from a set of initial parame- 
ters. Once the optimized parameters are determined, the 
expectation values are estimated with a large number of 
Monte Carlo steps. 

In the process of OWMC we must evaluated (cJ CT Cj CT ) to 
estimate the kinetic energy, which is the second process 
costing a lot of time in the simulation. We feel that there 
is a room to improve our algorithm on this point. 

Following a first step toward a development of the vari- 
ational theory based on the off-diagonal wave functions 
performed in this paper, we are planning to consider some 
problems for correlated electrons. In particular, a possi- 
bility of the superconductivity for the non-half-fillcd band 
is considered to be an interesting one. An application to 
other models also deserves an investigation. 

We thank Professor H. Aoki and Dr. K. Kuroki for 
useful comments and discussions. We express our sin- 
cere thanks to Professor P. Fulde for his suggestion about 
considering exp-5 type functions in the simulation. Our 
computations were supported by Research Information 
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